* Run event studies and DiD
* 2/9/2025

clear all
set more off

program main	
	* Static regressions:
	foreach design in sc state triple {
		twfe_regressions, design(`design') sample(actual_ssi)
	}
	static_csdid, design("sc") sample(actual_ssi)
	static_csdid, design("state") sample(actual_ssi)
	
	percent_effects
	generate_tables
	
	* Dynamic regressions:
	foreach design in sc state {
		dynamic_regressions, design(`design') sample(actual_ssi)
	}
end

program twfe_regressions
	syntax, design(str) sample(str)
	
	use "$for_analysis/`design'_for_regressions_`sample'", clear
	local indiv_controls = "i.sex i.race i.hispan i.edcat age i.diffphys"
	local state_snap_controls = "i.has_bbce i.has_call i.has_faceini i.has_facerec i.has_oapp i.has_no_fp"
	
	if "`design'" == "sc" {
		ren inter inter_sc
		lab var inter_sc "Single x Post"
		qui su snap if relyr == -1 & treat == 1, d
		global base_sc = `r(mean)'
		di "baseline SNAP takeup: `r(mean)'"
		reghdfe snap inter_sc [w = perwt], a(i.hhtype i.year) vce(cl statefip) // no controls
		est sto sc1
		qui estadd local cont " ", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype "\checkmark", replace
		qui estadd local stateFE " ", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_sc1 = _b[inter_sc] / $base_sc

		reghdfe snap inter_sc `indiv_controls' [w = perwt], ///
			a(i.hhtype i.year i.statefip) vce(cl statefip)	// controls
		est sto sc2
		qui estadd local cont "\checkmark", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype "\checkmark", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_sc2 = _b[inter_sc] / $base_sc

		reghdfe snap inter_sc `indiv_controls' ///
			state_unemp perc_democrat_2004 hh_income_2003 `state_snap_controls' ///
			[w = perwt], a(i.hhtype i.year) vce(cl statefip) // state controls instead of state FE
		est sto sc3
		qui estadd local cont "\checkmark", replace
		qui estadd local unemp "\checkmark", replace
		qui estadd local hhtype "\checkmark", replace
		qui estadd local stateFE " ", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_sc3 = _b[inter_sc] / $base_sc
	}
	else if "`design'" == "state" {
		lab var inter "Treat x Post"
		qui su snap if relyr == -1 & !mi(startyear), d
		global base_state = `r(mean)'
		di "baseline SNAP takeup: `r(mean)'"
		reghdfe snap inter [w = perwt], a(i.statefip i.year) vce(cl statefip)	// no controls
		est sto state1
		qui estadd local cont " ", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype " ", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_state1 = _b[inter] / $base_state

		reghdfe snap inter `indiv_controls' [w = perwt], ///
			a(i.statefip i.year) vce(cl statefip)	// controls
		est sto state2
		qui estadd local cont "\checkmark", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype " ", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_state2 = _b[inter] / $base_state

		reghdfe snap inter `indiv_controls' ///
			state_unemp `state_snap_controls' ///
			[w = perwt], a(i.statefip i.year) vce(cl statefip) // add state controls
		est sto state3
		qui estadd local cont "\checkmark", replace
		qui estadd local unemp "\checkmark", replace
		qui estadd local hhtype " ", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_state3 = _b[inter] / $base_state
	}
	else if "`design'" == "triple" {
		lab var triple "Single x Treated State x Post"
		qui su snap if relyr == -1 & treated == 1 & single == 1, d
		global base_triple = `r(mean)'
		di "baseline SNAP takeup: `r(mean)'"
		
		reghdfe snap triple [w = perwt], ///
			a(i.hhtype##i.year i.year##i.statefip i.statefip##i.hhtype) vce(cl statefip)	// no controls, only FEs
		est sto triple1
		qui estadd local cont " ", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype "\checkmark", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_triple1 = _b[triple] / $base_triple

		// note: state SNAP controls are collinear with i.year##i.statefip for triple diff, so need to omit
		reghdfe snap triple `indiv_controls' [w = perwt], ///
			a(i.hhtype##i.year i.year##i.statefip i.statefip##i.hhtype) vce(cl statefip)  // controls and FEs
		est sto triple2
		qui estadd local cont "\checkmark", replace
		qui estadd local unemp " ", replace
		qui estadd local hhtype "\checkmark", replace
		qui estadd local stateFE "\checkmark", replace
		qui estadd local yearFE "\checkmark", replace
		global beta_triple2 = _b[triple] / $base_triple
	}
end

program static_csdid
	syntax, design(str) sample(str)

	use "$for_analysis/`design'_for_regressions_`sample'", clear
	local indiv_controls = "i.sex i.race i.hispan i.edcat age i.diffphys"
	local state_snap_controls = "i.has_bbce i.has_call i.has_faceini i.has_facerec i.has_oapp i.has_no_fp"
	
	if "`design'" == "sc" {
		replace startyear = 0 if treat == 0 // csdid requires gvar = 0 for controls
	}
	else if "`design'" == "state" {
		replace startyear = 0 if startyear == .
		drop if statefip == 45 //drop South Carolina, as csdid excludes always-treated units anyway
	}
	
	* no controls
	csdid snap [w = perwt], gvar(startyear) time(year) notyet cluster(statefip)
	global N_`design'_csdid1 = e(N)
	estat simple, estore(`design'_csdid1)
	global beta_csdid_`design'1 = r(b)[1,1] / ${base_`design'}
	
	* with controls
	csdid snap `indiv_controls' [w = perwt], ///
		gvar(startyear) time(year) notyet cluster(statefip)
	global N_`design'_csdid2 = e(N)
	estat simple, estore(`design'_csdid2)
	global beta_csdid_`design'2 = r(b)[1,1] / ${base_`design'}
	
	* with state controls (cannot be time-varying)
	if "`design'" == "sc" {
		csdid snap `indiv_controls' perc_democrat_2004 hh_income_2003 [w = perwt], ///
			gvar(startyear) time(year) notyet cluster(statefip)
		global N_`design'_csdid3 = e(N)
		estat simple, estore(`design'_csdid3)
		global beta_csdid_`design'3 = r(b)[1,1] / ${base_`design'}
	}
end

program percent_effects
	di "baseline SNAP takeup rates:"
	di "sc: $base_sc"
	di "state: $base_state"
	di "triple: $base_triple"

	di "percent effects:"
	* twfe
	matrix percent_twfe = ${beta_state1}, ${beta_state2}, ${beta_state3}, ///
		${beta_sc1}, ${beta_sc2}, ${beta_sc3}, ${beta_triple1}, ${beta_triple2}
	mat li percent_twfe
	
	* csdid
	matrix percent_csdid = ${beta_csdid_state1}, ${beta_csdid_state2}, ., ///
		${beta_csdid_sc1}, ${beta_csdid_sc2}, ${beta_csdid_sc3}, . , ., .
	mat li percent_csdid
end

program generate_tables 
	* TWFE
	esttab state1 state2 state3 sc1 sc2 sc3 triple1 triple2, ///
		keep(inter inter_sc triple) order(inter inter_sc triple) ///
		varlabels(inter "Treated x Post" inter_sc "Single x Post" triple "Single x Treated x Post") ///
		se(%5.3f) star(* 0.10 ** 0.05 *** 0.01) stats(yearFE stateFE hhtype cont unemp N, ///
		labels("Year FE" "State FE" "Household FE" "Individual Controls" "State Controls" "N" ) ///
		fmt("%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc")) scalars() ///
		nomtitles label nonote booktabs
		
	esttab state1 state2 state3 sc1 sc2 sc3 triple1 triple2 ///
		using "$output/takeup_twfe", replace ///
		keep(inter inter_sc triple) order(inter inter_sc triple) ///
		varlabels(inter "Treated x Post" inter_sc "Single x Post" triple "Single x Treated x Post") ///
		se(%5.3f) star(* 0.10 ** 0.05 *** 0.01) stats(yearFE stateFE hhtype cont unemp N, ///
		labels("Year FE" "State FE" "Household FE" "Individual Controls" "State Controls" "N" ) ///
		fmt("%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc")) scalars() ///
		nomtitles label nonote booktabs
		
	* csdid
	esttab state_csdid1 state_csdid2 sc_csdid1 sc_csdid2 sc_csdid3 ///
		using "$output/takeup_csdid", replace ///
		keep(ATT) order(ATT)  se(%5.3f) star(* 0.10 ** 0.05 *** 0.01) ///
		stats(yearFE stateFE hhtype cont unemp N, ///
		labels("Year FE" "State FE" "Household FE" "Individual Controls" "State Controls" "N" ) ///
		fmt("%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc" "%9.0fc")) scalars() ///
		nomtitles label nonote booktabs
	foreach design in state_csdid1 state_csdid2 sc_csdid1 sc_csdid2 sc_csdid3 {
		di "for spec. `design', N = ${N_`design'}"
	}
end 

program dynamic_regressions
	syntax, design(str) sample(str)
	
	use "$for_analysis/`design'_for_regressions_`sample'", clear
	local indiv_controls = "i.sex i.race i.hispan i.edcat age i.diffphys"
	local state_snap_controls = "i.has_bbce i.has_call i.has_faceini i.has_facerec i.has_oapp i.has_no_fp"
	
	if "`design'" == "sc" {
		*TWFE:
		local omit = 6
		
		reghdfe snap ib`omit'.relyr_pos `indiv_controls' ///
			[w = perwt], a(i.hhtype i.year i.statefip) vce(cl statefip)
		est sto twfe
		
		*for csdid below:
		replace startyear = 0 if treat == 0
	}
	else if "`design'" == "state" {
		*TWFE:
		local omit = 6
		
		reghdfe snap ib`omit'.relyr_pos `indiv_controls' ///
			[w = perwt], a(i.statefip i.year) vce(cl statefip)
		est sto twfe
		
		*for csdid below:
		replace startyear = 0 if startyear == .
		drop if statefip == 45 //drop South Carolina, as csdid excludes always-treated units anyway
	}
	
	*CSDID:
	csdid snap `indiv_controls' [w = perwt], ///
		gvar(startyear) time(year) notyet cluster(statefip)
	
	estat event, window(-7 7) estore(cs) 
	
	coefplot (twfe, offset(-0.06) keep(*.relyr_pos) drop(1.relyr_pos 13.relyr_pos) ///
		rename(1.relyr_pos = Tm6 2.relyr_pos = Tm5 3.relyr_pos = Tm4 ///
		4.relyr_pos = Tm3 5.relyr_pos = Tm2 6.relyr_pos = Tm1 7.relyr_pos = Tp0 8.relyr_pos = Tp1 ///
		9.relyr_pos = Tp2 10.relyr_pos = Tp3 11.relyr_pos = Tp4 12.relyr_pos = Tp5) ///
		ciopts(recast(rcap) lcolor(blue)) mcolor(blue) ///
		mfcolor(white) msymbol(T) msize(medsmall) label(TWFE) baselevel) ///
		(cs, offset(0.06) keep(Tm5 Tm4 Tm3 Tm2 Tm1 Tp0 Tp1 Tp2 Tp3 Tp4 Tp5) ///
		ciopts(recast(rcap) lcolor(red)) label(DR) mcolor(red) mfcolor(white) msymbol(D) msize(medsmall)), ///
		xtitle("Years Since CAP Adopted") ytitle("Average Effect and 95% CI") vertical ///
		xscale(titlegap(2)) xline(5.5, lp(dash) lcolor(gs12)) yline(0, lwidth(thin) ///
		lp(dash) lcolor(black)) ylab(#15, angle(0)) graphregion(fcolor(white) lcolor(white) ///
		lwidth(vvvthin)) legend(region(lstyle(none))) ///
		coeflabel(Tm5 = "-5" Tm4 = "-4" Tm3 = "-3" Tm2 = "-2" Tm1 = "-1" Tp0 = "0" ///
		Tp1 = "1" Tp2 = "2" Tp3 = "3" Tp4 = "4" Tp5 = "5")
	graph export "$output/dynamic_0016_`design'.png", replace
end

* Execute
main
